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, Using equations of motion accurate to tiie third post-Newtonian (3PN) order (0(u/c)^ beyond 

CN| ' Newtonian gravity), we derive expressions for the total energy E and angular momentum J of the 

orbits of compact binary systems (black holes or neutron stars) for arbitrary orbital eccentricity. We 
j^', also incorporate finite-size contributions such as spin-orbit and spin-spin coupling, and rotational 

and tidal distortions, calculated to the lowest order of approximation, but we exclude the effects 
of gravitational radiation damping. We describe how these formulae may be used as an accurate 
diagnostic of the physical content of quasi-equilibrium configurations of compact binary systems of 
black holes and neutron stars generated using numerical relativity. As an example, we show that 
quasi-equilibrium configurations of corotating neutron stars recently reported by Miller et al. can be 
\ fit by our diagnostic to better than one per cent with a circular orbit and with physically reasonable 

^ . tidal coefficients. 
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I . I. INTRODUCTION AND SUMMARY 

m 

■ The late stage of inspiral of binary systems of neutron stars or black holes is of great current interest, both as a 
^ , challenge for numerical relativity, and as a possible source of gravitational waves detectable by laser interferometric 
O^' antennas. Because this stage, corresponding to the final few orbits and ultimate merger of the two objects into one, is 
I ^ highly dynamical and involves strong gravitational fields, it must be handled by numerical relativity, which attempts 
bJ[)' to solve the full Einstein equations on computers (see Refs. for reviews). 

^ ' The early stage of inspiral can be handled accurately using post-Newtonian techniques, which involve an expansion 
• i-H . of solutions of Einstein's equations in powers of e (^/c)^ ~ Gm/rc^, where v, m, and r are the typical velocity, mass 
' and separation in the system, respectively. By expanding to very high powers of e, one can derive increasingly accurate 
, formulae to describe both the orbital motion and the gravitational waveform. CurrentlVjresults for the orbital motion 
■ - - ■ accurate through 3.5PN order (0(e^/2) beyond Newtonian gravity) are known |3. II IH. M. II II liol ITll 111 111 IT^ . 

An important issue in understanding the full inspiral of compact binaries is how to connect the PN regime to the 
numerical regime. This is a non-trivial issue because the PN approximation gets worse the smaller the separation 
between the bodies. On the other hand, because of limited computational resources, numerical simulations cannot 
always be started with separations sufficiently large to overlap the PN regime where it is believed to be reliable. 
This has given rise to the so-called Intermediate Binary Black Hole (IBBH) problem for example, which seeks 
new techniques or insights to attempt to bridge the gap between the end of confidence in PN methods and the 
beginning of realistic numerical simulations. On the other hand, if it can be demonstrated that PN approximations 
converge sufficiently rapidly, especially for comparable-mass binary systems, then IBBH techniques may not be needed. 
Blanchet [ITgL [l7j | has recently argued that, for comparable-mass systems, the PN approximation seems to be more 
accurate than might be expected based on experience with the test-body limit. For binary neutron stars, this is less of 
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an issue, because neutron stars are much larger objects, so the numerical simulations necessarily commence at larger 
separations, where PN methods are presumably more reliable. 

Numerical simulations of compact binary inspiral start with a solution of the initial value equations of Einstein's 
theory; these provide the initial data for the evolution equations (some initial-data models solve in addition one 
of the six dynamical field equations). The initial state is assumed to consist of two compact objects (neutron stars or 
black holes) in an initially circular orbit. For stellar-mass systems that have evolved in isolation for eons, gravitational 
radiation is expected to leave the orbit in an accurately circular state, apart from the adiabatic inspiral induced by 
the loss of orbital energy; that inspiral is ignored in the initial-data models. (Miller has analysed the consequences of 
this particular assumption H^)- 

The circular-orbit condition is imposed by demanding that dr/dt = initially, where r is a measure of the orbital 
separation. One way to achieve this is to require that the system have an initial "helical Killing vector" (HKV), which 
corresponds to a kind of rigid rotation of the binary system. Some initial-data models assume that the objects are 
co-rotating, a condition which is astrophysically unlikely, albeit computationally advantageous, while others assume 
that the bodies are irrotational, i.e. non-rotating in an inertial frame. To simplify the problem, an approximation for 
the spatial metric is generally made; one is the assumption of conformal flatness, an approximation that is known to 
be invalid in full general relativity. This approximation is usually justified by the neglect of radiation reaction in the 
initial state. Other approximations, derived from post-Newtonian theory, or from sums of Kerr geometries, have also 
been used. For black hole binaries, suitable horizon boundary conditions must be imposed, while for neutron star 
binaries, equations of hydrostationary equilibrium and an equation of state must be provided. 

One important product of these initial value solutions is a relationship between the energy E and angular momentum 
J of the system as measured at infinity, and the orbital frequency fi. The energy could be the total energy as measured 
at infinity, consisting of the masses of the two stars plus the orbital energy, or it could be the total energy less the 
energy of the same two stars in isolation. The latter quantity would be a measure of the orbital binding energy. As all 
quantities are well-defined and gauge invariant, they are useful variables for making comparisons with PN methods. 

We have developed formulae for E{Q) and J(ri) using PN methods. Our analytic formulae include point-mass 
terms through 3PN order, but ignore radiation reaction. They also include rotational energy and spin-orbit and 
spin-spin terms for the case in which the bodies are rotating. They further include a Newtonian calculation of the 
effects of tidal and rotational distortions, applicable to stars of arbitrary density distribution, expressed in terms of 
so-called "apsidal constants" (i.e. we do not restrict attention to homogeneous ellipsoids HOI), and including effects 
at quadrupole and octupole order. We verify that, for black holes, tidal effects can be ignored, while for neutron-star 
binaries, they must be included. In contrast to previous work [T^ l2lll23. . our formulae apply to general eccentric 
orbits, not just to circular orbits. 

In an earlier paper , we compared this formula with HKV numerical solutions for corotating binary black holes 
obtained by Grandclement et al. for the regime where the black holes are separated from the location of the 

innermost circular orbit by a factor of around two, where PN results might be expected to work well (Gm/rc^ ^ O.I). 
We found that when we assumed circular PN orbits, our 3PN formulae for E{VL) and J{Vi) agreed to within 0.5 % with 
other PN methods, including our own formulae truncated at 2PN order, and 3PN formulae derived using resummation 
or Pade techniques. However all PN methods consistently and systematically underestimated the binding energy and 
overestimated the angular momentum, compared to the values derived from the numerical HKV initial-data models, 
by amounts that were up to 10 times larger than the spread among the PN methods. But when we relaxed the 
assumption of a circular orbit and demanded only that dr/dt — 0, our PN formula could be made to agree extremely 
well with the numerical data by assuming that the system being simulated is initially at the apocenter of a slightly 
eccentric orbit. For values of GmVl/c' ranging from 0.03 to 0.06, corresponding to orbital v/c between 0.3 and 0.4, 
or orbital separation between 10 and 6 Gm/c^, nearly perfect agreement with the binding energy and the angular 
momentum could be obtained with eccentricities that range from 0.03 to 0.05. 

The concordance within fractions of a percent between the various 2PN, 3PN and resummation PN results matches 
expectation, since (Gm/rc^ Y' ^ 10^"^. Presuming that all relevant physical effects have been included, we argued that 
the PN results in this range of GmVl/c' are robust. We suggested the possibility that the approximations made in 
most numerical initial-data models could lead to an apparent eccentricity in what was expected to be a quasi-circular 
orbit. At present, however, the discrepancy between the two approaches can only be considered a hint of possible 
eccentricity, because the results of Ts^ did not include quantitative error bars for the variables E{VL) and J{Vl). 

These results motivate us to propose a "post- Newtonian diagnostic" , a tool that can be used to extract physical 
information from numerical simulations, and that may also be an aid to guide some of the assumptions and approxi- 
mations inherent in numerical initial data computations toward those that lead to the desired physical configuration, 
such as a true quasi-circular orbit. 

In this paper we provide the physical assumptions, mathematical details, and justifications for the approximations 
that underly this proposed diagnostic tool. We give the detailed foundations for the analysis carried out in 24] for 
black-hole binary systems, and also extend that work to the case of neutron-star systems by including tidal effects. 
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As an application of our diagnostic to neutron-star systems, we analyse recent numerical models of quasi-equilibrium 
orbits of neutron stars by Miller et al. [251. In contrast to the black-hole case, we find that the orbital energy in 
the neutron-star initial-data models of can be fit to better than one percent, and importantly, within the error 
bars provided in using circular orb" with physically reasonable tidal parameters appropriate to the 'T = 2" 
equation of state used in that numerical work. The results illustrate the robustness of the PN approximation well into 
the strongly relativistic regime of compact binaries, especially when augmented with physically movitatcd finite-size 
effects. Application of this PN diagnostic to other numerical models will be subject of future papers. 

The remainder of this paper provides the details underlying these conclusions. In Sec. ^1 we solve the post- 
Newtonian equations of motion calculated to third post-Newtonian (3PN) order, for general eccentric orbits. Neglect- 
ing radiation reaction effects, we then express the total conserved orbital energy and angular momentum in terms of 
a pair of "covariant" orbit elements e (eccentricity) and C, (related to the semi-latus rectum). In Sec. IIIII we calculate 
the effects of finite size in binary systems with bodies whose spin axes are perpendicular to the orbital plane. These 
include tidal and rotational distortions, spin-orbit terms and spin-spin terms. In Sec. IIVI we analyse our diagnostic 
quantitatively, and apply it to co-rotating, equal-mass binaries of black holes and of neutron stars. Two appendices 
provide the detailed derivations of the expressions for the tidal and rotational distortion included in our diagnostic: 
Appendix ^ uses Newtonian gravity to solve the general problem of the equilibrium configurations of gravitating 
bodies disturbed by an external force, paralleling the treatment in the classic monographs of Kopal [2a. |2^. and 
Appendix ^ specializes the results to linear perturbations caused by rotational and tidal disturbances. 



II. ENERGY AND ANGULAR MOMENTUM FOR "POINT" MASSES TO 3PN ORDER 



A. Orbits at the turning point in post-Newtonian gravity 



Since our ultimate focus will be on orbits that are possibly eccentric, but that momentarily have r = 0, it will be 
useful to review the characteristics of orbits at turning points in Newtonian theory. In Newtonian gravity, the orbit 
of a pair of point masses may be described by the set of equations 



p/r 

E 
J 



1 + e cos(0 — Lo) , 
r^d(l)/dt = [mpY'^ , 

^|x X v| , 



(2.1) 



where p = a{\ — e^) is the semi-latus rectum (a is the semi-major axis), lo is the angle of pericenter, m = mi + m2 is 
the total mass, ^ — 11111112 /m is the reduced mass, and E and J are the total orbital energy and angular momentum, 
respectively (henceforth we use units in which G = c — 1). A circular orbit corresponds to e = 0, with r = a = 
constant, Vl^ = m/a^, E/fi = a^f2^/2 — m/a ~ — (mf2)^/^/2, and J//i — ^Jrna = {m/il)^^^. However, if we demand 
only that the orbit be at apocenter, so that r — only, we have 
that, in terms of fla, the angular velocity at apocenter, 



LU + TT, r = p/{\ — e), f]-^ ~ {m/p^){l — e) 



J / fim 



(l-e)2 

1-1/3 



2/3 



(1 



(2.2) 



flp and 



To obtain expressions in terms of fip, the angular velocity at pericenter, one makes the replacements fla 
e — e in Eqs. H2.2|l . 

However, at higher PN orders, neither the orbital eccentricity e nor the semi-latus rectum p is uniquely or invariantly 
defined. One definition of eccentricity used by Lincoln and Will in their analysis of orbits at 2.5PN order was 
that of a Newtonian orbit momentarily tangent to the true orbit (the "osculating" eccentricity); it had the unusual 
property that it did not tend to zero for a circular PN orbit, but tended toward a constant value of order m/p, while 
the rate of pericenter advance approached the same rate of rotation as the orbit itself. In this language, the true orbit 
was a non-circular orbit at perp etual periastron, thereby maintaining a constant separation r. In an effort to avoid 
this anomaly, other authors [Sfll adopted a "quasi-Keplerian" parametrization, which defined multiple "eccentricities" 
to encapsulate different aspects of non-circular orbits at PN order. 



4 



In an effort to find a parametrization of non-circular PN orbits that will be useful in comparing with numerical 
models, we p3 | proposed an alternative measure of eccentricity and semi-latus rectum according to: 



4/3 



C-- - ^ j , (2.3) 

where flp is the value of O where it passes through a local maximum (pericenter) , and Qa is the value of Q where it 
passes through the next local minimum (apocenter). 

These definitions have the following virtues: (1) they reduce precisely to the normal eccentricity e and semi-latus 
rectum p in the Newtonian limit, as can be verified from Eqs. 1)2. l|l : (2) they are constant in the absence of radiation 
reaction; (3) they are somewhat more directly connected to measurable quantities, since fl is the angular velocity as 
seen from infinity (eg. as measured in the gravitational-wave signal) and one calculates only maximum and minimum 
values, without concern for the coordinate location in the orbit; and (4) they are straightforward to calculate in a 
numerical model of orbits without resorting to complicated definitions of "distance" between bodies. 

They have the defect that, when radiation reaction is included, they are not local, continuously evolving variables, 
but rather are some kind of orbit-averaged quantities (for this reason, they may not be as "covariant" as they seem 
- see Sec. Ill El below). Nevertheless, when an eccentric orbit decays and circularizes under radiation reaction the 
definition of e has the virtue that it tends naturally to zero when the orbital frequency turns from ocillatory behavior 
to monotonically increasing behavior (i.e. the maxima and minima merge). 

By virtue of these definitions, C has the further property that 

O \ 2/3 / o \ 2/3 



We will derive expressions for orbital energy and angular momentum in terms of these parameters e and ^; for 
comparision with numerical models of quasi-equilibrium parametrized in terms of 17 at r = {fla or Op), one can 
simply substitute for ^ from Eq. 1)2.4(1 . In this section we will focus on 3PN expressions for point masses; in the next 
section, we will incorporate effects due to rotation and finite size. 



B. 3PN equations of motion 

We use the standard form of the equations of motion, written in a "Newtonian-like" manner. The acceleration of 
body 1 is given schematically by 

ai = ^^ = ^ {n[-l + (PiV) + (2PiV) + (2.5PiV) 

+ {3PN) + (3.5PiV) + ■■■] 
+v[{PN) + {2PN) + {2.5PN) 

+ {3PN) + {3.5PN) + ■■■]} , (2.5) 

where Xq and uia denote the position and the mass of the body a, r is the separation between the two bodies, 
n = (xi — X2)/r is the unit vector from 2 to 1, and v = vi — V2 the relative velocity. The equation for body 2 is 
obtained by making the replacement 1^2. The notation nPN represents the n*'' post-Newtonian correction to 
Newtonian gravity. These equations are valid only for point-like, non-spinning bodies. 

Post-Newtonian terms nPN include even (integer) and odd (half-odd integer, such as 2.5PN, or 5/2 PN) orders. 
Even terms are conservative, in the sense that the equations of motion admit conserved quantities such as energy and 
angular momentum. Odd terms correspond to gravitational radiation reaction, and therefore are not conservative. In 
particular, they will cause the orbit to shrink, and the eccentricity to decrease. 

We convert the two-body problem to an effective one-body problem. For this purpose we choose the origin to be at 
the center of mass of the system, which is defined by an integral of the motion (a conserved quantity to the 3PN order 
of approximation to which we will be working) . We then change all variables to the relative coordinates x = xi — X2 
using relations of the type 



Xi = [m2/m+ {riSm/2m){v'^ - m,/r) + {2PN) ~\ ]x, 

X2 = [—mi/m + {ri6m/2m){v'^ — m/r) + {2PN) + ■ ■ ■]x , 



(2.6) 
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where rj = /x/m = mim2/m? is the reduced mass parameter {0 < rj < 1/4), and Sm = mi — m2- The result is a set 
of equations of motion in terms of relative coordinates: 



(2.7) 



where A and B represent post-Newtonian terms. To date, the two-body equations of motion have been computed up 
to and including 3.5PN order. In an appropriate harmonic gauge, writing A = Ai+ A2 + ■ ■ ■ and B = Bi-\- B2 + ■ ■ ■ ^ 
the expressions for A and B read [3ll |: 
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(2.9a) 
(2.9b) 

(2.9c) 



(2.9d) 



(2.9e) 



At 3PN order, the computation implemented by Blanchet et al. j^, |9(| produced logarithmic terms, proportional to 
ln(r/r'i) and ln(r/r'2), where r'l and r'2 are constants related to a scale of radius for each body. In obtaining Eqs. 



(12. 8|) and (|2.9|) . we removed these logarithms using a 3PN coordinate transformation 



Sx^, with [^: 
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Sxf_, = — —mim2df. 
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where ija = |x — Xq| denotes the coordinate separation between the considered point and the body a. We note that we 
have □(Sx^ ~ 0, except at the location of the two bodies. This ensures that the harmonic condition is still respected 
in the new gauge to the required order. In addition, the parameter A, which was initially undetermined in |^ 0, |^ 
has now been fixed to be A = —1987/3080 by different techniques [Tllll^ls^: that value has been incorporated into 
all equations. 

In the absence of the 2.5PN and 3.5PN terms, these equations of motion admit conserved total energy E and total 
angular momentum J. Writing E = Eq -\- Ei -\- E2 + Ey, and J = Jq + Ji + J2 + J3, we have: 
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C. Solution of the 3PN equations of motion 



In order to solve these equations, we shall initially adopt the method of osculating orbital elements, which is well- 
adapted to the perturbed two-body Kepler problem. The osculating orbit elements are defined by the Keplerian 
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orbit that is tangent to the actual trajectory at a particular moment of time. In the Newtonian case, the oscula ting 
elements are constants of the motion; in a perturbed Newtonian problem, they change smoothly with time (see |29| | 
for more details about the method of osculating elements applied to the post- Newtonian problem) . 

From the equations of motion we can easily deduce that the trajectory is planar, which allows us to reduce the 
number of variables from six to four. If we assume that the plane of the motion is perpendicular to z (x,y,z being 
a standard cartesian coordinate system), our new set of variables {a,P,p,(j)) is related to the old set {x,y^Vx,Vy) by 
the definitions (some of which are redundant): 



X = r cos ( 



y = r sm 4> , 

Vx = -(rn/p)^/^(/3 + sin(?!.) , 

Vy = (m/p)^/'^(a + COS0) , 

r = p(l + a cos (/) + /3sin c/))^"'" , 

r^<P = (mp)i/2. (2.13) 

Reciprocally, we can deduce the osculating elements from the orbital variables by using the following relations: 
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One additional expression will be useful: 



(2.14) 



(to/p) (a sin 0-/3 cose/)). (2.15) 



We note that the vector (a, (3) has as its norm the ordinary Keplerian osculating eccentricity e and as its phase angle 
the direction u of the Keplerian osculating periastron, so that we have a = ecosw and /3 = esinw. 

In what follows, we will use the parameter u — m/p rather than p. Note that u is of order e ~ m/r. In the 
Newtonian case, u, a and (3 are constants of the motion; in the post-Newtonian problem, these parameters vary 
according to the following "Lagrange planetary equations" (so-called from their extensive use in solar-system studies) : 

— = -2^3/2^, 
dcp 

^ = As\n<l) + 2u^/'^B{a + cos4>), 
d(p 

^ = -Acos(j} + 2u^''^B([3 + si\Y(j)), (2.16) 

where we have used Eqs. (|2.7|l . H2.13|l and 12.15|1 . When the definitions of x and v [Eqs. (|2.13|) ] are substituted into 
the PN expressions for A and B [Eqs. H2.8I) and (|2.9|l ]. we get a set of coupled first-order differential equations in the 
variables a{4>), f3{4') and u{(j)). 

The planetary equations derived from Eqs. H2.15|) are too long to be reproduced here (they can be found through 
2.5PN order in '2^). However we can schematically write them in the general form: 

^ = u^Vui{a,P,(j)) +u^Vu2{a,P,(j}) +u^/'^Vu5/2{a,P,(j}) -\ , 

dq) 

— uVaila, + u^'Da2(a, f3,(t)) + u^^'^'Da5/2(a, (3,(1)) + ■ ■ ■ , 

dcp 

^ = (a, /?,(/.)+u2l?/32 (a, /3,0) +7/^/215/35/2 (a, /3, 0) + ••• , (2.17) 
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where Pun, 2?a„ and 'DPn (n G {1, 2, 5/2, 3, 7/2}) are polynomials in a and /3 and simple trigonometric functions of 
(j). We quote, for illustration, the first post-Newtonian expressions for these polynomials: 

Vui = 4{2 — T]) {/3 cos (j) — a sin (j>) , 

Vai = -3/3+ (3 - 77) sin./) + (5 - 477) (a sin 20- /3 cos 20) 

+ i[(56 - 47?7)a^ - [% + 2lri)0^] sin0- 7(32- 13r;)a/3cos0 
8 4 

3 3 
+ -ri{0^ - a^) sin 30 + -?7a/3 cos 30 , 
8 4 

= 3a- (3 -77) cos 0- (5- 4r;) (a cos 20 + /? sin 20) 

-i[(56 - 47?7)/32 - (8 + 21'q)a^] cos0 + ^(32 - 1377)0/? sin 
8 4 

3 3 
+ o'?("^ - ^^)cos30+ -77a/3sin30. (2.18) 
8 4 

We want to solve these equations iteratively. At zeroth (Newtonian) order u, a and /? are constants of the motion 
M, a and /3, and can be related to the initial state of the orbit. Post-Newtonian effects cause them to vary slowly 
over a post-Newtonian timescale or a radiation-reaction timescale, related to the orbital phase by 0/e and 0/e^/^, 
respectively. Superimposed upon this will be variations on an orbital timescale. To take these two effects into account, 
we use a two-scale approach j33j . We define a variable 9 — e0, and we assume that the osculating elements can be 
written as functions of 9 and in the generic form u = u{u{9), a{9), I3{9), 0), with 9 and now treated as independent 
variables. We then expand the elements in powers of e: 

u = + e^ui (a, /3, u, 0) + 6^^2(5, /3, u, 0) + ••• , 

a = a + eai{a, j3,u,(f>) + e^a2{a, (3, u,(f>) + ■■■ . (2-19) 
Notice that, by its very nature, u begins at order e. We write the derivative with respect to in the form 

d d 



(2.20) 




da d df3 d du d \ 

d9da ^ de'dp d9m J ■ 

We also expand the derivatives with respect to 9 in powers of e: 

du " ~ 

— — dui(a, P,u) + edu2(a, f3,u) + ■ ■ ■ , 
d9 

^ = d/3i(a,/3,-u) + ed/32(a,/3, u) -I , 

dct ~ ~ 

— = dai{a, P,u) + eda2ia, I3,u) + ■ ■ ■ . (2-21) 
du 

Now we have reduced our study to the search for at, Pi, Ui on the one hand, which will give the dependence on 
a, [3, u, and 0, and dai, d(3i, dui on the other hand, which will give differential equations allowing solution for the 
^-dependence, or long-term variation of the parameters. Note that this is not the only way to decompose the problem, 
but is a natural way, given the split into orbital and secular evolution of the variables. 
We now define the average and the average- free part of a function /(0) by: 



<.27r 
JO 



if) ^ i- 1 fm^, 



AT{m) = f[4')-{f). (2.22) 

where the "independent" variable 9 is held fixed. (An equivalent procedure would be to convert all functions of 
into 27r-periodic functions and constants.) We rewrite Eqs. (|2.1t)|) with our new variables, and we collect terms of 
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common powers of e. At first order in e we get 



dui + -^-j- = u^'Dui{a, (3, (/)) , 



dai 



dai 



uDai(a, (3, 4>) 



d/3i + ^ = uD(3i(6i,i3,^ 



(2.23) 



where the expressions on the right-hand-side are given by Eqs. H2.18|l . with a replacing a, and so on. Reading off the 
average parts of Eqs. H2.23|l . we find dui = 0, dai = — 3{t/3, and d(3i = Swa. Defining a = ecosw and $ = esinw we 
find, to first PN order that 



du/de = 0, 

di/dO = {&d&/d9 + Pd(3/d9)/i = 0, 
dOj/dO = {ad(3/de - f3d&/de)/i^ ^3u. 



(2.24) 



These results express the well-known fact that the orbital eccentricity and semi-latus rectum do not evolve secularly 
to fPN order; in fact, this holds true at 2PN and 3PN order; they only evolve secularly as a result of radiation 
reaction. The angle of pericenter Cj evolves secularly at f PN order via the standard advance; there are also 2PN and 
3PN contributions, but no radiation-reaction contributions to the advance of a), through 3.5PN order. 
Then, integrating the average- free parts of Eqs. (|2.23|) . we obtain, for example, 

(2.25) 



ai = AJ'{Vai){<j>)d(f) 



The role of the second AJ- is to get rid of the constant of integration. The same method yields similar results for jSi 
and ui. 

At second order in e, we obtain equations of the form 



dai 9ai ~ dai 

+uiDai — dai —dpi — dui 

da dp ou 

f2{a,P,u, (j)), 



(2.26) 



where ai, /3i, ui, dai, d(3i and dui are known from the first order solution. For the same reasons as previously we 
have: 



da2 = (/2) 
a2 = Af{ I AF(/2)(0)d0) 



(2.27) 



Using this procedure systematically up to 3.5PN order, we completely determine a{a, (3,u,(j)), P{a, (3,u, (/)) and 
M(a,^,'u, 0), as well as ^(a,/3, u), ^{a,f3,u) and ^{a,P,u). From this and Eqs. 12.1311 we can deduce the ex- 
plicit expressions for x, v, r, etc. 

To 3.5PN order, the secular evolution of u and e is governed by radiation reaction, and is given by the coupled 
equations (we now set e = 1) 



du 



de 



= iv\ (8 + 7e2){t7/2 



/2759 



V 42 



6?7 



/379 127 \ o A483 79 \ . 



-—T]e\ (304+ 121g2)u5/2 _ 
15 t 



/ 18049 



63677 



/4346 1829 



V 7 



(2.28) 
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We note that the eccentricity decreases as the orbit shrinks. The periastron advance is driven by the conservative 
part of the equations: 

duj 3 r , , ,T , fST /157 123 .\ , 

— = 3u--[l0 + 477-(l + 1077)g2]fi2^|__ (^___^2j^_3r?2 



45- 23 



123 
128' 



93 2 

TT" ] 1] + —f] 



+ -77(12 - 25ri)e'^ } 



(2.29) 



D. Energy and angular momentum in terms of new orbit elements 

We now wish to convert from the osculating orbit elements u and e to our alternative quantities defined in Eqs. 
(cf. Seeing. Using the formula 



mfl = m^/^p^/^/r^ = u^/^(l + a cos (/) + /3 sin ( 



(2.30) 



we can easily show that the maxima and minima of Q occur at (f> — uj (pericenter) and (p ^ lj + tt, (apocenter) 
respectively. We then express ftp and fla, and thence our new orbit elements e and ^ as functions of e and u. To 2PN 
order, the relationships are given by 



e~<!l 
1 



i(13-4r;) + (l-3r7)g2 



i(52 - 129ry) 



— (157 - 337?7 + 116ri^)e^ + -(4 - 19r/ + A8ri^)i^ 



C = u U - 3 [(3 - ^) + (1 - 3?7)e'] u + 



-(198 + 3977 + 26772) 



- — (1092 - 97777 + 27677^)^2 + 1(2 - 2777 - 18r]^)i^ 
36 9 



(2.31a) 



(2.31b) 



Notice that a circular orbit corresponds to e = e = 0. 

We invert these relations and substitute the expressions for e(e,C) and {((e,C) into the solution of the equations of 
motion. The results for m/r and r'^cj) to 3PN order are too long to be reproduced here. However, in order to give an 
idea of what they look like, we quote them to 2PN order, expressed in terms of our new orbit elements. 



y = |l + ecos0jc+|l-^'7+^(4-37;)e2 



+ - [(9 - 477) + (1 - 377)6^] ecos(/)'- |e2cos(20')|c^ 

1 - ^77 + ^17(356 - 31977 + 4877^)62 + -1:7(256 - 26577 + 45977^)6'' 

4:(96 - 2317/ + 877^) + -^(323 - 35I77 + 180rj^)e^ + ^{5 + 127j)e'^ 
12 48 12 



(60 + 1597/- I6772 



48 



(29 - 27r])e' 



cos(2(/)') 



1 

"I6 

TO 



(1 + 1977 - 477^)6^ cos(30') - ^(1 - 377)6^ cos(4</>') } ( 
1 - 1^(3 - 77) + ^(1 - 377)e2 + (4 - 27/)ecos(/.' \ C 



+ <{ 1(6 + 537/ + 277^) - ^(28 + II777 - 12772)6^ + 1(2 - I777 + 67]^)e^ 



1 



1(6 - 5377 - 277^) - ^(32 - 2117/ + 547/2)e2 



+ ^(36 - 137/ + 477^)6^ cos(2^') - ^(3 + 277)6^ cos(30') \ C 



ecos( 



(2.32) 
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where (j)' = (j) — uj. The leading term corresponds to the Newtonian solution. Note that C: e and cu are now our 
post-Newtonian orbital elements, and should not be mistaken for the Newtonian u, e and uj introduced in Sec. Ill Al 
An alternative method for integrating the post-Newtonian equations of motion was developed by Wagoner and Will 
|28|. In that method, perturbations of the velocity and angular momentum were defined by the equations 



r^d^/dt = \yiy.^r\ = {mpY''^{l + 5h) , 

V = (m/p)^/^[— sin((/) — w)ex + (e + cos(0 — aj))ey + ^v) 



(2.33) 



where p, e and w are constants. Taking a time derivative of both equations, substituting the 3PN equations of 
motion (ignoring radiation reaction terms), converting to derivatives with respect to cj) and integrating, one obtains 
expressions for the perturbed r'^d(p/dt and v. One then integrates the identity {d/d(j))r^^ — —{r'^dip/dt)~^{n ■ v) with 
respect to (j), setting the constants of integration at each PN order so that the identity d{rn)/dt = v is reproduced. 
In terms of the bare orbit elements e, u = m/p and lu, the orbit equations look different at each PN order from those 
derived above in terms of e, u and a). But when the solution so derived is used to identify fla and Qp and thence to 
define new orbit elements from Eqs. (|2.3(l . the resulting orbit solution in terms of our new orbit elements is identical 
to Eqs. through 3PN order. 

The equations describing the evolution with time of our new orbit elements then become: 



dC 



de 



du 



1186 685 



743 

h 22?7 

42 ' 



63 



6 



/ 18001 163 
V 1008 ~ ~6" 



^9/2 



-^r;e|(304+ 121e2)C^/2 
/ 12499 17741 



5505 3796 



3C 



V 21 
1 
2 

137 



6 

(9 - 1477) 
/337 



-V 
1 

^4 
123 
128 



7 

/ 46289 
V 168 



(19- 18r;)e2 
53 



3 

437?/^ e* 
27 

y ~ 



481 



123 



77 + 777 



TT^ I 77 



-'7 



e^ + -(20 + 8r7 + 4577")enC 



(2.34) 



Now the problem is entirely solved. Equations (|2.32|l pushed to 3PN order, characterize the motion, while Eqs. 
give the pericenter advance and effect of radiation reaction on the orbital elements. 

We now ignore the effects of radiation reaction, express all the orbital variables r, x, v, r to 3PN order in terms of 
our new orbit elements and the angle (j), and substitute into the expressions H2.11|l and H2.12|l . As expected E and J 
are constant (independent of </>) through 3PN order. Defining E = E/fi and J — | J|//i777, with J ~ | J|z, we find for 
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a general eccentric orbit: 
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24 

853 
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77 

432 ' 

1 
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5 

16 

„3 



3013 

7 
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(2.35a) 



(2.35b) 



Notice that i?Harm is proportional to (1 — e^) through 3PN order, indicating that SHarm ~ for the limiting unbound 
orbit e — 1; this is another appropriate feature of our "covariant" eccentricity. The energy and angular momentum are 
well-defined, physically observable quantities, so one can alternatively express our orbit elements C and e as functions 
of E and J. Here we give the results to IPN order, but the calculation can be done to 3PN order: 



c 



1 

J2 



1 



2 
3J2 



= Vl + 2E.P 



4 + 277 -(l-3?;)£^j2 
E 



1 

1 - - 



2 1 + 2EJ^ 



4 + 277 - (1 - 377)£;j^ 



(2.36) 



E. ADM vs. Harmonic gauge 

The foregoing results are valid in harmonic gauge. That gauge is characterized by the condition d^h^'^ = 0, where 
h'^'^ = ^/—gg^'^ — 77'''', g^ii/ and g are the physical metric and its determinant, and 77^^ is a background Minkowski 
metric. In this gauge, Einstein's equations take the form 

□V = le^T^,., (2.37) 

where □ — 77'^''i9^9j^ is the fiat d'Alembertian operator, and the source term t^^, depends both on the matter stress- 
energy tensor T^^ and on non-linear contributions of the gravitational field. The local equation of motion Wi.T'^" = 
is equivalent to d^r'^'^ = 0, which follows from the harmonic gauge condition. There actually is an infinity of distinct 
harmonic gauges, and the equations of motion will generally depend on the choice of a particular gauge. We already 
saw an example of this in the choice of eliminating logarithmic terms from the 3PN contributions [Eq. (|2.10|l above] . 

A difi'erent approach to the two-body problem, implemented through 3PN order by Damour, Jaranowski and Schafer 
0,0,15, is to compute the Hamiltonian of the system rather than the equations of motion. Unlike other methods, 
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3r), 

m 



this does not use a harmonic coordinate system, but a so-called ADM ( Arnowitt-Deser-Misner) , or "Hamiltonian" 
gauge, or coordinate system. It has been proven to be equivalent to the harmonic formulation 

The Hamiltonian has been computed up to 3PN order; because it is a Hamiltionian approach, it explicitly suppresses 
the 2.5PN and 3.5PN contributions of radiation reaction. We quote it here only to IPN order (see for a complete 
expression) : 



Ha 



DM 



Pi 

2toi 



mim2 
2r 



_j1 

8ml 



miTO2 / (n • pi)(n • P2) 



4mim2 



2m\ 



7 Pi ■ P2 

4 77117712 



-(1 



We convert this two-body problem into an effective one-body problem by using the simple relation p 
valid in the center-of-mass frame. Thus we get a new expression for Habm (x, p) . 
From Hamilton's equations: 



Pi 



dx 
It 



Vpi?ADM 



dp 



ADM, 



(2.38) 



P2, 



(2.39) 



we iteratively extract the equations of motion and write them in the same form as equation H2.7|l . but with different A 
and B. Substituting the expression of p as a function of v and x into the Hamiltonian, we obtain the total conserved 
energy i^ADM- Similarly, we get Jadm by calculating x x p. For both the equations of motion and the expressions for 
energy and angular momentum, the harmonic and ADM-Hamiltonian terms coincide at IPN order, but they differ at 

2PN and 3PN orders. 

We apply the method described in section III CI to find solutions to the ADM equations of motion and expressions 
for -Badm and Jadm in terms of the osculating orbit elements. In this case, e and u are strictly constant because 
radiation reaction is not present in the Hamiltonian approach. We then find expressions for our new orbit elements e 
and C in terms of e and u and write -Eadm and Jadm in terms of these elements. The results are: 
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We observe two features of the harmonic and the ADM versions of these expressions: (i) the "circular" parts (e — 0) 
of the formulae coincide. In that case the angular velocity = ila = ^Ip is the same as that observed from infinity 
for both harmonic and ADM coordinates; (ii) the expressions also coincide for 77 ^ 0, i.e. in the test-mass limit. As 
mentioned before, the differences between the formulae only occur at 2PN and 3PN orders. It is actually possible to 
relate the coordinate positions and velocities in the two gauges. In particular, the relation between (/jadm and (jjuarm, 
?'Harm, etc. allows US to find a relation between (badm, Cadm) and (eHaim, CHaim), and thus account for the differences 
in the coefficients of E and J. We found that a transformation of the type 



0ADM = 




(2.41) 



where we have dropped the subscript "Harm" in the right-hand-side of Eq. (|2.41|) . and where / is a function, was 
compatible with the differences observed in the expressions of both the energy and the angular momentum. Since 
f = at the apastron and periastron, / does not need to be determined explicitly for our purposes. In the circular orbit 
limit, where, from Eq. (|2.7|) . = m/r[(l — (3 — ri)m/r] to PN order, it is easy to see that 0adm — '/'Harm- Eq- (|2.41|) 
demonstrates that our definitions of e and C are not truly covariant. Nevertheless, the coordinate transformations that 
connect different formulations of the post-Newtonian equations of motion cause changes beginning only at 2PN order. 
This is reflected in Eq. H2.41|) where the difference between the two angular velocities is of 2PN order. Furthermore, 
for the small eccentricity orbits that we wish to consider, the corrections are proportional to e, and are thus further 
suppressed. Thus we argue that our definitions of e and C are "almost" covariant. 



III. EFFECTS OF FINITE SIZE 
A. Estimates for compact binaries 



In reality, the bodies in our binary system cannot be treated as purely point masses. They may be rotating, and thus 
subject to a number of effects, including rotational kinetic energy, rotational flattening, and spin-orbit and spin-spin 
interactions. Furthermore, there will be tidal deformations. These effects will not only make direct contributions to 
the energy and angular momentum of the system, they may also modify the equations of motion, and thereby modify 
the expressions for our alternative eccentricity and semi-latus rectum. However because they depend on the size of 
the bodies, which, for neutron stars and black holes, are of order m, we expect these effects to be "effectively" of high 
PN order, even if they are Newtonian in origin, such as tidal effects. To see this, we estimate each finite-size effect in 
turn and compare it with the Newtonian orbital energy En ^ jr. We assume that the rotational angular velocity 
Lo of each body ranges from zero to the orbital angular velocity, given by 17 ~ {m/r^Y^'^ , and we let the radius of each 
body be of the form Ra ^ quia, where q 1 for black holes (in harmonic coordinates), and g 5 for neutron stars. 

• Rotational kinetic energy: i?R,ot ^ Iuj"^ /2 < mR^{m/r'^) ~ Efqq^im/rY . This is effectively 2PN order. There 
will be PN corrections to the kinetic energy, given by i?Rot-PN ~ -BRot(-R^)^ ^ mR'^uj'^ ^ E]yq'^{m/r)^ . These 
are effectively 5PN order, but, because of the g"* dependence, could be important for neutron stars. 

• Rotational flattening: iSpiat ^ 6Ilu'^/2, where 5 is a measure of the deformation of the body, given by the ratio 
of rotational to gravitational energy, S ~ {luj^) / {m? / R), so that iJpiat ~ ijo'^R^ < Eis[q^{m/r)^ . There is an 
equivalent contribution of rotational flattening to the gravitational internal energy. These are effectively 5PN 
order, but because of the q^ dependence, could be important for neutron stars. 

• Tidal deformations: E'Tidai ~ {5'm)'^/R, where 6' is the ratio of gravitational energy due to the tidal force 
of the companion to the internal gravitational energy of the body, 5' ~ {mR^ / r^) / {m / R) ~ (R/r)^. Thus 
^'Tidai ^ m?{R/rY' / R ^ ENq^{m/r)^ . There is also a contribution from the rotational kinetic energy of the 
tidal bulge, given by i?KE-buigc S'Iuj"^ ~ mR^uj'^/r^ < E^q^ {m/r)^ . These are effectively 5PN order, but 
could be significant for neutron stars. 

• Spin-orbit couphng: -Es.o. LS/r^ ^ {mr'^VL){mR^uj)/r^ < ENq^{m/rY . This is effectively 3PN order [3^. 
and generally must be included. 

• Spin-spin coupling: E's.S. ^ 5'i S2 Ir^ ~ {mR^u>Y /r^ < EMq^im/r)^ . This is effectively 5PN order, but could be 
significant for neutron stars |34| . 
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A parallel heirarchy of finite-size effects applies to the total angular momentum of the system. 

The largest effect in principle is that due to the rotational kinetic energy of the bodies and thus requires some care. 
For black holes, we can apply the general formulas for mass and angular momentum of isolated Kerr black holes, in 
terms of the irreducible mass and angular velocity. For neutron stars, no such general formula exists, so it may be 
necessary to rely upon numerical results for energy and angular momentum of isolated rotating neutron star models 
in order to take accurate account of this effect. On the other hand, it does not directly affect the equations of motion. 

Because the remaining effects are effectively of 3PN order and higher, our strategy will be to evaluate them analyti- 
cally to the lowest non-trivial order. For tidal and rotational flattening terms, this will mean using Newtonian theory. 
For spin-orbit and spin-spin terms, we will use the well-known IPN formulae. We will ignore any coupling among 
these effects, or between these effects and the point-mass PN effects described in the previous section. Accordingly, 
we will calculate the separate contribution of each effect to the energy and angular momentum and simply add them 
all up. 



B. Newtonian Tidal and Rotational Effects 



In Appendix ^ we derived the general form of the equations of motion and the conserved energy and angular 
momentum for a binary system of tidally and rotationally deformed bodies, and in Appendix IBl we specialized to 
linear perturbations and multipole indices I = 2 and 1 = 3. We now specialize further to systems more relevant to the 
initial configurations in numerical relativity which we wish to study, namely binary systems in which the spin axes of 
both stars are perpendicular to the orbital plane. The equation of motion (|B9d|l then takes the simplified form 



1 



r / \ r 



(3.1) 



where the three perturbing terms correspond repectively to the effects of rotational distortions, quadrupole tidal 
distortions {I = 2) and octupole tidal distortions {I = 3), with the coefficients given by 



A 
B 
C 



8m^'^ (Rlk^^^ 1712 / mi + R2k^^^mi/'m2) ■ 



(3.2) 



For each body, Ra denotes its radius, k2'^^ and fcg"'' denote the "apsidal constants" for angular harmonics I = 2 and 
I = 3, respectively, and (Da denotes the body's angular velocity at a chosen point in the orbit (see Appendix IbI for 
details). Apsidal constants are dimensionless coefficients that depend on the degree of central condensation of the 
star, and that determine the size of distortion of a given angular degree I produced by a given external perturbation. 



Note that A < R^/m^ 



q°{m/ry, so that, despite appearances, this term, like the purely tidal term from B, is 



effectively 5PN order. The energy and angular momentum that are conserved by virtue of the full fluid equations of 
motion are given by 
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J- 



TR, Orbit 



[IlUJl + I2CU2] 

X v| , 



(1 



(3.3) 



where, for each body, la denotes the moment of inertia, Wa denotes the self-gravitational energy of the undistorted 
configuration, and f denotes the orbital separation at the point at which the star's angular velocity is Co. The chosen 
point in our case will be the pericenter or apocenter. In Eq. (|3.3|l . the split among the intrinsic energy and spins of 
the bodies i^seif and S, the constant distortion terms -^Distort and JDistort, and the orbital terms is clear. The angular 
momentum components are all referred to the axis perpendicular to the orbital plane. 
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We now repeat the method of subsections III CI - III Dl to obtain the general solution to the equations of motion to 
first order in the tidal and rotational perturbations. We then obtain our new orbit elements e and C in terms of the 
bare elements e and it; for example, e is given by 



e = e<l-- 



1 + ^e~2 



8 8 



1 

~ 4 
931, 
256^ 



, 85,2 85,4 
1 e e 



12 



Ci 



24 



(3.4) 



Since we are assuming that these effects are effectively of 5PN order, we can simply add the correction terms in Eq. 
(|3.4|) to those in Eq. (|2.31a|) . Tidal and rotational interactions are conservative (as long as we ignore dissipative 
processes such as viscosity), and therefore do not cause secular evolution of e or u; however they do produce a 
pericenter advance, given in terms of our new orbit elements by 



— ^ AC^ 



. 3 2 1 4 



BC 



15 



15 



1 + — + —e^ + —e 
4 8 64 



(3.5) 



Substituting the solutions for the motion into the orbital parts of Eqs. 1)3. 3|l . and converting to our new elements, 
we obtain for the tidal-rotational (TR) contributions to the orbital parts of E and J, 



E' 



TR, Orbit 



i(3 - e2)^C^ + J-(9 + lOe^ - 3e^)BC'^ 
9 18 



+ :t7(13 + 49e2 + 7e'^ - 5e^)CC* 
24 



'^TR.Oi 



bit 



^(3 + 6^)^(3/2 + ^(3 + lOe^ + 3e4)BC'/2 



+ -(l + 7e2 + 7e'^ + e^)CC^^/^ 



(3.6a) 



(3.6b) 



where we have dropped the Newtonian orbital part, because it is already included in the 3PN point-mass expressions 
of Eqs. (|2.35|l or (|2.4U|) . The form of the self-terms depends on where in the orbit we evaluate the stars' angular 
velocities; for pericenter or apocenter, we can use the Newtonian relation that m/f = ({1 ± e), respectively, to write 

- Wi + (1 ^ 2) , (3.7a) 
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^-Distort 
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imuJi 



)2 + 2^C'(l±e)^l +(1^2) 



2{mCjif + 3— C^(l ± ef] +(1^2), 
m i 



(3.7b) 
(3.7c) 

(3.7d) 



C. Spin-orbit and Spin-Spin effects 



Spin-orbit and spin-spin interactions produce corrections in the equations of motion that are formally of IPN order. 
For systems with the spins perpendicular to the orbital plane they are given by 



1 



6n(n X v) • 2S 



V X 7S 



6m . 
3 — A 
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(^3S4 








TO J } 



rnSi • S2 



(3.8) 



where S = Si + S2 and A = m(S2/TO2 — Si/toi). The individual spins are constants of the motion when they are 
both aligned perpendicular to the orbital plane. The conserved energy and total angular momentum are given by 



1 2 1 T 

E = -^v h -^Ln 



J = L 



N 



s-v 
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S + — A) - lsi.S2, 
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(3.9a) 
(3.9b) 
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where Ljv = iix x v, and Eq. Ij3.9b|l denotes the component perpendicular to the orbital plane (for the complete 
equations of motion, see, for example, [ssl l36l|l. We define the dimensionless quantities 



D = ri- 



S 



F = ri- 



5m A 
m L 



N 



G = r] 



S1S2 



(3.10) 



where D ^ F ^ (R/r)'^ ~ q'^{m/r)'^, and G ~ {R/r)^ ~ q^{m/r)'^, making the spin-orbit and spin-spin terms 
effectively 3PN and 5PN order respectively [s^. With these definitions, the equation of motion takes the form of Eq. 
(E3, with 



A = {5D + 3F - 3G)v^ ~ 3{D + F - G)r^ 
B = -2Dr. 



(3.11) 
(3.12) 



Again we solve the equations of motion using the method of subsections III CI - III DI and define our new orbit elements. 
In this case, for example, the eccentricity is given by 



= g |l i [(1 -f 4i^)D + (3 + 2e^){F - G)] uj 



In terms of our new elements, the pericenter advance is given by 

^ = -(7£i + 3F-3G)C. 



(3.13) 



(3.14) 



while e and ( undergo no secular changes. When expressed in terms of our new orbit elements, the spin-orbit and 
spin-spin contributions to the total energy and angular momentum have the form 



1 



- e") [(7 - 26^)0 + (3 - ^)[F - G)] , 



Jspin = -iM"^[5(7 + e2)i?+(15 + e2)F-4(3 + e2)G] VC- 
Inserting the Newtonian expression for Ljv, we have that 

D^{Slm^)^, F ^ i6m/m){A/m^)^, G ^ iSi/mi)iS2/m2)m-\ . 



(3.15a) 
(3.15b) 

(3.16) 



D. Other finite-size corrections 



In deriving the "point-mass" equations of motion, the underlying assumption was that the masses that enter the 
equations are the total mass of each body, comprised of baryonic mass, gravitational binding energy and rotational 
kinetic energy, if any. Thus, each nia should be written iria = rn^ — Wa + E^°^. In many numerical approaches, 
sequences of models are constructed in which the total (or ADM) mass of each corresponding non-rotating star 
is held fixed along the sequence. Thus, for making comparisons with such sequences, we should replace each ma 



in Eqs. or with 



laio\l2 (or, in the case of black holes, with a suitable formula in terms of 



the irreducible mass and lOo). But because la^\ ~ q^imlrf"^ the main contribution, at effectively 3PN order, 
comes from making this replacement in the Newtonian expressions. Expressing Ejsi and J^r in terms of f2 as Ejsi — 
— ^rim{l — e^){mQa)^^'^ / ~ e)*^'^, and J^v = r\rr?{^ — e)^!"^ j {m£La)^^'^ ^ and making the above replacement, we find 
the corrections to the Newtonian energy and angular momentum 



^^N,Corr 



fim 



mi 



-?(l-fi)+(1^2) 
1 V 3m/ 

(1^2) 



nil \ 
3m J 



(3.17a) 
(3.17b) 



where all masses now are those of the equivalent non-rotating body. For neutron stars, this would be that of the 
same baryonic mass; for black holes, it would be that of the same irreducible mass. 
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Equal Mass vs. Extreme Mass Ratio 

Contributions relative to Newtonian orbit energy 




FIG. 1: Contributions of IPN, 2PN and 3PN terms to the energy, expressed as a fraction of the Newtonian energy, vs. mfl. 
Circular orbits are assumed. Shown are the equal mass case {rj — 1/4) and the point-mass limit [i] — 0). 



IV. A POST-NEWTONIAN DIAGNOSTIC FOR QUASI-EQUILIBRIUM CONFIGURATIONS 

A. Estimates of effects 

We now have all the ingredients to formulate a post-Newtonian diagnostic for quasi-equilibrium configurations of 
compact bodies. The ingredients are the various contributions to the total energy and angular momentum of the system 
in terms of the "covariant" orbit elements e and together with the relationships connecting the value of C. with the 
orbital angular velocity at a turning point of the orbit, namely C — (mr2p)^/^/(l + e)''/'^ or — (mr2a)^/'^/(l — e)^/^, 
corresponding to pericenter and apocenter, respectively. The ingredients are: 

• Point-mass orbital contributions through 3PN order. Eqs. (|2.35|) or (|2.4U|I . It is straightforward to show that, 
because the harmonic and ADM versions differ by 2PN terms proportional to rje^ and higher, the differences 
between the two versions are negligible for all cases of interest. Henceforth we will adopt the harmonic version 
of Eqs. EHSJl. 

• Self Terms. Eqs. H3.7a|l and (j3.7b|l . We add a suitably defined total "rest" mass for the bodies to the definition 
of £'scif- Because the rotational kinetic energy and the spin angular momentum are effectively of 2PN order, 
they will have to be treated with some care. 

• Constant distortion terms. Eqs. (|3.7c|) and Ij3.7d|l . 

• Tidal-rotational orbit terms. Eqs. H3.6|l . 

• Spin- orbit and Spin- spin terms. Eqs. H3.15|l . 

• Newtonian correction terms. Eqs. H3.17|l . 

In order to assess the applicability of this diagnostic, we first study the sizes of various effects for systems of interest. 
In general we will consider systems of solar-mass scale neutron stars or black holes, in circular or small-eccentricity 
orbits, in the vicinity of the onset of an unstable plunge and merger. This corresponds to C~m/r< 1/5 for black 
holes, or to C < rn/{2R) ^ 1/q for neutron stars. For q between 4 and 6, the two ranges are comparable. Both 
correspond to toJ7 < 0.1. We will generally choose a range 0.01 < mfl < 0.1. 
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Corotating Neutron Stars (R=4m) 

Contributions relative to Newtonian orbit energy 
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FIG. 2; Contributions of tidal and spin terms to the energy for corotating neutron star binaries, expressed as a fraction of the 
Newtonian energy, vs. mf2, for q = Ra/ma = 4. Circular orbits are assumed. Shown are the PN contributions for comparison 

First we look at the relative contributions of point-mass PN corrections. Figure 1 shows the contribution, relative 
to the Newtonian orbital energy, of the IPN, 2PN and 3PN terms in the energy, for 77 = and 77 = 1/4, as a function 
of mf7. Results for the angular momentum are similar. While the IPN terms are essentially insensitive to ry, and the 
2PN terms are only 15 % smaller for equal masses than for the test-mass limit, the 3PN terms are suppressed for equal 
masses by more than a factor of 10 compared to the test-mass limit. As Blanchet 0,^3 has argued, this suggests 
that the 3PN approximation may be quite accurate for comparable-mass systems, without the need for sophisticated 
resummation techniques. At the largest angular velocity considered, 3PN terms contribute less than one per cent of 
the total binding energy and angular momentum of the orbit. 

Next we consider the effects of tidal and rotational distortions. We consider systems of identical bodies (mi = m^) 
which are corotating {Cji — Cj2 — For neutron stars, we adopt the maximum values of the apsidal constants 
{k2 = 3/4 and = 3/8, see Appendix IB 3|l . and choose two representative values of g = Ra/ma for neutron star 
models with reasonable equations of state, namely g = 4 and q = 6. The results, plotted as a fraction of the Newtonian 
orbital terms, are shown in Figures 2 and 3, along with the PN contributions for comparison. As expected, tidal effects 
are very sensitive to the stellar radii. For g = 4, the I = 2 tidal terms become comparable to the 2PN and IPN terms 
only around mQ ~ 0.09, while the I = 3 terms are an order of magnitude smaller. For q = 6, the I = 2 tidal terms 
exceed the IPN terms already by mfi ~ 0.05, while the I — 3 terms are small, approaching the 2PN terms only at the 
largest allowed mil ^ 0.07, corresponding to the point at which these larger stars are touching. For irrotational stars 
(wi = (I'2 = 0), the tidal effects are very similar. 

These curves illustrate that tidal effects need to be taken into account carefully in an accurate diagnostic for neutron 
star binaries, but are not so large that they invalidate our approximation scheme. Their modest size also supports 
our use of Newtonian theory to calculate them. They only become problematical for the largest neutron stars near 
the very endpoint of their inspiral. It should also be pointed out that, in making these estimates, we have adopted 
the largest values of the apsidal constants, corresponding to uniform-density stars. While neutron stars are not as 
centrally condensed as, say, non-degenerate stars, they are also not uniform density, so the fc; may well be smaller 
than their maximum values. For example, for a Newtonian polytrope, p = kp^ , with F = 2, A:2 = 0.26, so the q = 6 
tidal terms in Fig. 13 are reduced by a factor of three, bringing them to a level at or below the IPN terms over the 
whole range of mfl. On the other hand, very little, if anything, is known about the values of fc; for general relativistic 
neutron stars over a range of equations of state. This is a subject that we are currently investigating. 

Figure 4 shows the effects of tides for co-rotating black-hole binaries. There we choose q — 1 (i? = to in harmonic 
coordinates), k2 =3/4 and k^ = 3/8 (for slowly rotating black holes, fc2 from rotational distortions happens to be 
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Corotating Neutron Stars (R=6m) 

Contributions relative to Newtonian orbit energy 

\ I I r 




FIG. 3: Same as Fig. ^ with q = 6 




FIG. 4: Same as Fig|21but for corotating black-hole binaries, with q — 1. 



precisely 3/4; see, eg. We see, not surprisingly, that tidal effects are utterly negligible over the entire range of 

mfl. 

Finally, we examine spin effects. Again we consider identical, co-rotating bodies. For neutron stars, we assume that 
Sa — la^, with the moment of inertia given by that for a uniform density body, la — {2/5)maR^ = {2/5)q^ml^. The 
results are shown also in Figures 2 and 3. For g = 4, spin-orbit effects are small but significant, just below the 2PN 
terms, while spin-spin effects are negligible. For q = 6, spin-orbit terms exceed 2PN terms by mfl ~ 0.04 and become 
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comparable to the IPN terms by the maximum angular velocity, while spin-spin terms barely exceed the 3PN effects. 

For black holes, we use the fact that Sa — 4m^r2. Figure 4 shows that the spin-orbit terms lie between the 2PN 
and 3PN contributions and thus must be included, while spin-spin terms are negligible (though larger than the tidal 
terms). 

B. Corotating, identical black holes 

For black hole binaries, we ignore tidal and spin-spin effects. We set mi = m2, rj = 1/4, and uji — 1I12 = ^■ 
We exploit the fact that there exist exact formulae for the energy and spin of isolated Kerr black holes in terms of 
the irreducible mass, M = Mi„/[1 - 4(Mtow)2]i/2^ s = iM^„uj/[l - 4(Mirrw)^]^/2. The total energy and angular 
momentum of the system are then given by 



^Tot — -Eself + i?Harm + £'N,Corr + -Espin , 
JtoI — S + Jnarm + >/n,Coit + >/Spin , 



(4.1) 
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(4.2a) 
(4.2b) 

(4.2c) 

(4.2d) 
(4.2e) 
(4.2f) 
(4.2g) 
(4.2h) 



where mi„ is the total irreducible mass of the system, given by (mii-i)! + (miii.)2. In Eqs. (|4.2a|l and Ij4.2bp . we have 
expanded the Kerr formulae for M and S in powers of mi,-,-fl, assumed to be small compared to unity, keeping as 
many higher-order terms as needed to reach a precision comparable to our 3PN formulae. To obtain Etoi and Jxot at 
a turning point as functions of fi, we substitute C — (mirrria)^^"^/ (1 — e)^^'^ or ( = (mji-rilp)^/'^/ (1 + e)"'/'^ for apocenter 
or pericenter, respectively (in calculating i?N,Corr and JN.Com we have already changed the dependence in C from 
the total mass of the rotating bodies to the total irreducible mass of the non-rotating counterparts). These are the 
formulas used in [23| to compare with the numerical HKV quasi-equilibrium solutions of Grandclement et al. [l^ . 
When i?Tot and Jxot are scaled by mi„ and m?j, respectively, there remains only one free parameter, the eccentricity 
of the orbit, and we found .24;] that a substantially better fit to the numerical data was obtained for non-zero values 
of e, of the order of 0.03, with the system at apocenter, than for e = 0. We suggested that such apparent eccentricity 
could be a result of the inevitable approximations (such as the conformally flat approximation) and numerical errors 
in such initial-data models, but, in the absence of detailed estimates of the sizes of those errors, it was difficult to 
draw firm conclusions. On the other hand, those engaged in numerical models of black hole binaries could use our 
diagnostic as a guide to know when, say, a suitable circular orbit has been achieved, or whether further numerical 
experiments with different grid sizes or larger computational domains are necessary to reach the desired physically 
meaningful state. 
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C. Corotating, identical neutron stars 



For neutron stars, we must include tidal effects. We set mi — TO2, 77 = 1/4, and cDi 0)2 = 51; we let the apsidal 
constants and radius factors be common for both stars, given by k^, and q, respectively, and express all quantities 
in terms of the total mass tuq = (mo)i + (mo)2 of two non-rotating stars with the same equation of state. We also 
define for each star the coefhcient aa = la/maR^, and also assume it to be common for both stars. The result is 
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where the 3PN point- mass expressions i?Harm and Jnarm are given in Eqs. H4.2c|l and (|4.2d|l . and where 
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We illustrate the use of this diagnostic by comparing with numerical data recently reported by Miller et al. |25| . 
They constructed a sequence of general relativistic, quasi-equilibrium configurations of corotating neutron stars, in 
the conformally flat approximation. They used a polytropic equation of state with F = 2. Among other quantities, 
they report an "effective" binding energy, given by E), = [A/adm — '2'M^q,{^)]IMq, as a function of Afofi, where Madm 
is the total ADM mass of the configuration, and MNs(fl) is the ADM mass of a uniformly rotating isolated neutron 
star of the same baryonic mass Mq, as each star in the binary configuration, but rotating with angular velocity 57. 
Since the rotational kinetic energy of the stars is already removed, we can compare the numerical results with the 
FN diagnostic i?Diag = [ETot — -Escifl/TOQ. Since our toq is twice the ADM mass of a non-rotating neutron star, 
we must scale E^, by Mo/2AfADM-NSi where A^adm-ns is the ADM mass of an isolated, non-rotating neutron star. 
In the models of Miller et al., AfADM-NS = Afo/1-067. We also need to fix the coefficients q and a. From data 
provided by Miller (private communication), the radius of each isolated non-rotating star in isotropic coordinates is 
given by Ri = 6.77A/ADM,while the baryonic moment of inertia, calculated using isotropic coordinates, is given by 
/o = 9.412A/o. We work in harmonic coordinates, but since Rh = -R/(l + AfADM/^^/); the difference between the 
two coordinates is only of order 1/2 %, so we read off q = 6.77. The ADM moment of inertia can be identified as 
/adm = (AfADM/A/o)/o = 9.412AfADMA'/o^ = 9.412(1.067)2A'/|j3m- Thus we can read off aq^ = 9.412(1.067)^ and 
hence a — 0.234, or around half of the uniform-density value of 2/5. (Miller also calculates the same quantities in 
terms of circumferential, or Schwarzschild radius; after transforming to harmonic coordinates, the results for q and a 
are consistent with these to within a few per cent.) 

Inserting these values of q and a into our diagnostic, we compare various FN configurations with those reported by 
Miller et al., as shown in Figure 5. The numerical results are shown as with error bars, estimated by Miller et 
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Binding Energy of Corotating Neutron Stars - Numerical vs. PN 
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FIG. 5: Comparison of PN diagnostic with numerical initial-data models of |2^. 



al. from the results of a range of convergence tests. Four curves show the energy for circular orbits, for various values 
of the apsidal constants. Neither the uniform-density values {k2 = 3/4, ks = 3/8), nor the point mass values (fe — 0, 
^3 = 0) gives a good fit at all, except at low angular velocities (large separations) where tidal effects are smaller, 
and all circular-orbit curves converge toward the numerical result. Models with half the uniform-density values for 
k2 and k^ give marginal fits. However, a very good fit is achieved with values fc2 = 0.260 and k^ = 0.106; these 
are precisely the values for Newtonian F = 2 polytropes (Appendix , which is the equation of state used in the 
Miller et al. numerical models. Also shown is a model with the same F = 2 apsidal constants, but with a non-zero 
eccentricity e — 0.02 and with the system at apocenter. This marginally fits the numerical data within the error bars, 
but consistently gives lower (more negative) energies. 

We conclude that these quasi-equilibrium neutron-star configurations are fit to better than one per cent by our PN 
diagnostic with a circular orbit, and with physically reasonable tidal terms. 

In future work we plan to compare this diagnostic with results of other numerical models of quasi-equilibrium black 
hole and neutron star binaries. Our 3PN equations of motion, together with tidal and spin terms, augmented by 
radiation reaction terms, can also be used to develop a "d yna mical" diagnostic, to compare with numerical simulations 
of evolutions from the quasi-equilibrium initial data pH l38| | . 
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To derive the effects of tidal and rotational flattening, we will adopt standard methods from Newtonian theory for 
binary systems, such as those detailed by Kopal |2^ l27j. We assume that the timescale for changes in perturbing 
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APPENDIX A: NEWTONIAN TIDAL AND ROTATIONAL EFFECTS 



1. Distorted equilibrium configurations 
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quantities (such as the external tidal potential, seen either from the global inertial frame, or from the rotating 
frame of a given body) is sufficiently long that each body can be assumed to be in hydrostatic equilibrium. In 
other words, we will ignore dynamical tides j39l |. This is a reasonable assumption as long as we are focusing on 
quasi-equilibrium initial data. Consider one of the bodies in the binary system. From the equation of hydrostatic 
equilibrium, Vp — pV'I', where p, p and ^' are the pressure, density and total gravitational potential, respectively, we 
conclude that Vp x V^* — 0, and thus that surfaces of constant p and ^E" coincide. We label surfaces of constant p by 
the radial parameter a, and let the equation of those surfaces have the form 



ria. 



e,ci>)=a[l + Y,fiMYi„,{(l)], (Al) 



where Yim{^) are spherical harmonics corresponding to the direction fi, and where the dimensionless distortion 
functions have the property /*^ = 

On general grounds we expect fim ^ {R/rY^^ ^ (7'+^(m/r)'+^ for tidal effects, and, for I = 2, /2m ^ t^^/p ^ 
{R/r)^ ^ q^ijn/rY from rotational effects. The effect of these distortions on the external potential of a body is of 
order fimiR/fY ~ 5^'"'"^(m/r)^'+^. For I — 2, this means effectively 5PN order; I — 3 effects would be effectively 
7PN order, and so on. However, for neutron stars, with (7^4 and m/r ^ 0.1, an ^ = 2 distortion effect becomes 
numerically comparable to a 2PN term, while Z = 3 is comparable to a 3PN term. For black holes, with q < 1, the 
effects are much smaller. Thus, in the end, we will keep only I = 2 and I = 3 distortion terms. Also, non linear 
corrections to fim would be of order (R/rY^^ ^ q'-^^{m/rY^^ smaller than the dominant linear effects, and thus, 
effectively of 8PN order for I = 2 (for neutron stars, these non- linear corrections would be numerically smaller than 
3PN) . The exception to this is in the internal gravitational energy of each body, where a quadratic contribution yields 
(jn^/R)fim^ ~ {m? /r){R/rY^^'^^ , which is comparable to the other effectively 5PN contributions for I — 2. 

We begin, however, with a general analysis, keeping Im arbitrary, and working to second order in the small quan- 
tities /;,„. Later (Appendix IB|) we will speciaHze to I = 2 and I = 3 hnear perturbations. To second order, it is 
straightforward to show that, for any n. 



r" = a"{l + n}^[/,™(a) + (n-l)X,™]rz„(r!)}, (A2) 

where 

= 2 51 ^'^■nsfafifjS , (A3) 
and Cl^.^g is defined in terms of Clebsch-Gordan coefficients, 
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Note that the various angular momentum quantum numbers are connected by the constraints I — a + "/ , a + 
7 — 2,... ,\a — 7I, and m = l3 + S; the C^^.^g are symmetric under {af3) ^ (jS). Also note that Xqq — 
(16^)-V2^^^/„^/*^. 

We expand the gravitational potential U of the body and the disturbing potential V in the form 



Air r r 

^ = y p(x')^r,:„(J7<)yz™(f7>)dV, 



A-TT 

V = dr^ + J2——di^r'Yi^{n) , (A5) 

Im 

where the subscript > (<) corresponds to the larger (smaller) of r and r' . The disturbing potential consists of a part, 
with disturbing coefficients dim, that corresponds to a potential with V^y — 0, such as the gravitational potential from 
another body, or the Laplacian-free part of a centrifugal potential, plus the spherical part of a centrifugal potential, 
with coefficient d. We now substitute Eqs. IjAip and ljA2p into (|A5|I . convert all expressions from r to a, and demand 
that, for I = 0, the external gravitational potential of our body have the form U = m/r (i.e. the perturbation does 
not change the mass of the body), and that, for I 7^ 0, the total potential U he constant at a given a. The first 
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can be satisfied if /oo + 2Xoo = 0, while the second holds if, for I ^ 0, 



Att a 

El 21 + 1 



-fim{a) + a dim 



2^ + 1 m{a) 
27r a 



Xim{a) ~ 2{2l + l)da^ fi^ , 



(A6) 



where 



m{a) 
Fim{a) 
Eim{a) 



At: p[a)a? da , 



Z"' d 



p{a)da 



da 



a'+^fim + {l + 2)Xim,) 



(A7) 



and A denotes the surface of the body. The left-hand-side of Eq. ljA6|l is first order in while the right- hand-side 
is second-order. Dividing the first-order terms by a', differentiating with respect to a and multiplying by a^'+^, we 
obtain the first-order result 



(A8) 



where prime denotes differentiation with respect to a. Substituting this and the first-order solution of Eq. (jA6l) 
back into the right-hand-side of Eq. (|A6(I . it is straightforward to show that the term involving Yl'afS-yS reduces to 
m{a){2Xim — aX[^) / A'Ka, to second order. The basic equation for the distortion functions can then be written in 
the form of an integral equation 



{a)fim{a) - a-' / p{a'+^ fim)' da - a'+^ / p{a^-'f,^yda 



21 + 1 
An 

a'+M;„, + 7^^,„(a) , 
where TZim contains all contributions quadratic in small quantities: 

21 + 1 



(A9) 



Att 



am(a)XL + 2(2^ + l)daV/^ 



+ {l + 2)a-' / pia^+^XhnYda 



-(/-l)a'+i / p{a^-'Xir„yda. 



(AlO) 



Combining Eq. ljA9|l with various derivatives of it, one obtains the following useful equations, evaluated at the 
surface a = A of the star: 



Att 

il + l)firniA)-AfUA)+Vira = — A"' / p[a'+^ {fim + (l + 2)Xim)]' da , 



lfl,n{A) + AfLiA) - Qir. 



An 



(Alia) 
(Allb) 



where m — Anpa^da, and 



8n ~ 



Vim = A^X\'„,{A)~lAXUA) + —dAMfUA)~{l~l)fUA)], 

m 

Qim = A^Xl'^{A) + {l + l)AX[^{A)] + —dA^[AfUA) + {l + 2)fim{A)] 

m 



(A12) 
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Another combination of first and second derivatives of Eq. (jA9|) yields a second-order differential equation for 
sometimes called Clairaut's equation: 

a'/L + ^^{afL + firn) - l{l + 1)/;™ = -^yn'L - l{l + 1)%™] . (A13) 
mya) m\a) 

For a given density distribution p{a)^ this equation can be solved, subject to the boundary conditions that Jim be 
regular at a = 0, and that, at the surface, fim satisfy Eq. (|Allb|l . 



2. Energy and angular momentum of the system 

Given a solution for the distortion functions fimia), we can calculate all the quantities needed for the equations of 
motion and the energy and angular momentum of the orbit. The external potential of our body, for example, is given 
by 

A 



U = ^+E'^^/ p[a'-^(/,™ + (/ + 2)X.„)]'da 



"m sr^' 47r Yim(^)mA , , , , , ,^ 



r 

I 



7 + E ^i:^^JTrf'^rndlmYlm{n) , (A14) 
Im 



where denotes summation for / ^ 0, and where we define the "apsidal constant" ki^ by 

l{l + l)fi^(A)-Afl^{A)+Vim 



2 lfi„XA)+Af[^^{A)~Q,rn 
The total gravitational energy of the system is given by 



The self-energy Wn of body 1 can be written 



Substituting Eqs. HA1|I and (|A2|) we find no contribution linear in To second order, we obtain 

Jo a 

\a^\fL\^ + 2afUL + + 



(A15) 



W = -\j j ^-^^d'xd'x' = W,, + 14^12 + (1^2). (A16) 



^11 = -E^^ / pr^^'drdnYimm F p'r''+^dr'dn'Y*^in') . (A17) 

;™ ^« + t Jq Jq 



(A18) 



Since the second term is already second order, we can integrate by parts and use the first-order versions of Eqs. (jAllj) 
and IjAlSp to obtain the alternative form 

W,,^-W + J2'^^A^'^'\dim\'kim, (A19) 

Im 

where we define the self-gravitational binding energy W of the undistorted configuration by 

>V= r'^dmia). (A20) 
Jo a 
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In Wi2, we substitute the external potential of body 2 evaluated inside body 1, to obtain 

(2) (2) 

where y2 = x — X2, is the apsidal constant of body 2, and is the coefficient of the disturbing potential acting 
on body 2. Since the interaction energy is smaller than the self energy by a factor of R/r, we only need to keep terms 
linear in the deformations fim or the disturbing coefficients dim', consequently we carry out a multipole expansion of 
1 /?/2 in the first term in Eq. l)A2Hl , then convert from r to a using Eq. (EH, but we evaluate the second term at the 
center of mass of body 1 and do the lowest-order spherical integral. Effectively, we are ignoring multipole-multipole 
coupling between the bodies, which can be shown to lead to effects of order (A/r)^'+^(A/r)^"+^ ~ (m/r)^°, or lOPN 
order for Z = n = 2. The result is 

Im 

where now r = |xi — X2I and n = (xi — X2)/r. Combining Eqs. ljA19|) and (|A22p . the final result is 

1 mi 7719 

2 r 
47r 



Im 



+ (1^2), (A23) 

where, under the interchange, n — > — n. 

The kinetic energy of the system is given by T — / pv^d^x. Splitting the velocity of an element of fluid into 
center-of-mass, rotational, and random parts, and noting that 

E YiMYUn') = = ^^'(" • ' (A24) 

m 

where n"^^^ denotes an STF product of I unit vectors (a capitalized superscript denotes a multi-index), the product 
n^^^n'^^^ denotes contraction on all indices, and Pi is a Legendre polynomial, we may write 

T = + T^l,^^, + \loI P pr'drdn ( 1 - y E Y2,n{Xi)Y,*^m] +(1^2), (A25) 

where Ai = iVi/ui. Converting from r to a using Eq. IjAljl . recalling that /oo — — 2Xoo, and noting that is already 
of first order in disturbing quantities, we obtain, to second order in small quantities, 

T = Im^vf + r«_, + \hul - ^^u^lAl ^ d'ilk^^lY.miXi) + (1^2), (A26) 



where Ii = (2/3) Airpa^da. 

The angular momentum of the system is given by J* = e^^^ j px^v^. Using the same split of the velocities, we 
obtain, to the analogous order of precision, 

,r = mi (xi X vi)^ + hio\ - 2iu{A\ S^lk'ilz<^> + (1^2), (A27) 

where we define the symmetric trace free (STF) tensor 

Z<^>= fYim{n)n<'^>d^n, (A28) 



''Ini 

with the following properties 



A^^,<f> = j2TTTju^'"'^^^^- ^^^^^ 
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3. Equations of motion 



The Newtonian equations of motion for body 1 are given by 

a\ = — [ pdPx I p'V- — ^—d^x'. (A30) 
mi Ji J2 |x-x'| 

We write x = xi + x and x' = X2 + x' and expand in a Taylor series about xi and X2. We define the STF muhipole 
moments 1^'°^^ — px'^^^d^x with /° = TOq and = 0. Finally, we calculate the relative acceleration a* = a\ — a\. 
After some manipulation, we obtain the general result 

mr* ^1 f T<^> T<^>\ /r 



a = — h m 



Y.n[- — + (-1 ^ 



^ ^ Z! \ mi 7712 / \r 



+"^> > 1/, M — ^ - ' (A31) 

where m = mi + 1112 and the products of the multipolc moment tensors are to be symmetrized on all indices and 
made trace- free. For our distorted bodies, the STF multipole moments can be shown to be given by 

/<^>=2Af+i^d«e)z<f>. (A32) 

m 

With the coefficients dim ^ m/r'-^^, we have that /<^> ~ TOy4^'+^/r'+^, and therefore the multipole-multipole 
coupling term in the equation of motion ljA31|) is of order (TO/r^)(A/r)^'^+^; since (? > 4, this is lOPN and higher. As 
before, we ignore multipole-multipole terms. 



4. Multiple disturbance sources 



We will want to consider both tidal disturbances as well as rotation-induced disturbances. To see how this affects 
our general results, we note that the non-linear corrections to the Clairaut equations never play a role to the order of 
accuracy we require, only the linear functions //,„, satisfying linear differential equations, are needed in the end. Let 
fim — 9im + him , whcrc each disturbance function satisfies the linearized Clairaut equations ljA13p , with a boundary 
condition for each determined by the linearized Eq. (|Allb|l . From the structure of the formulae for the external 
potential U, the kinetic energy, the angular momentum, and the multipole moments, it is clear that the coefficient 
dimkim can simply be replaced by d^'^k^^ + d\'^k'l^, where is the amplitude of the disturbing function for that 

disturbance, and is the corresponding apsidal constant, determined from the linearized Eq. (|A15p . However 
because it has a contribution quadratic in disturbing functions, the gravitational self-energy Wn requires some care. 
Returning to the expression (jAlSp . substituting fi^ = gim + ^/m and carrying out the integrations by parts, using the 
linearized Clairaut equations satisfied separately by gim and him, one can show that the coefficient \dim\'^kim must be 
replaced by the coeflicient \4^2 1'^l^J + \d^J I'k^'J + kl^Mt^^J:^ + ^.^^dl^^^i^). 



APPENDIX B: ROTATIONAL AND 1 = 2, I = Z TIDAL DISTORTIONS 



1. Disturbing coefficients and apsidal constants 



We focus on the lowest-order I = 2 and I = 3 tidal terms. The gravitational potential at a point x' in body 1 due 
to body 2 is given by 

, ^ 4_ 

V = ^2Y^ ^f:^^i-^yYl*rnin)Ylm{n') , (Bl) 
Irn 

where we ignore the contributions to V from the distortions of body 2 (ignore multipole-multipole coupling). We thus 
obtain the tidal coefficient for body 1, 

dj^'^ =m2i~iyYCmin)/r'+\ (B2) 
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with the coefficient for body 2 obtained by interchange. Working to hnearized order, we can factor out the azimuthal 

fT _ fTy* I 
llm — Jl ^lm\ 



m dependence by defining — ffYimin), then the Clairaut equation and the outer boundary condition for body 1 



take the form 



a'ff" + + fF) l{l + l)fl = , (B3a) 

A,friA,) + lfnA,) = 4^-1/^ (^y^' . (B3b) 

Note that the apsidal constant depends only on ff, and is thus independent of to, so 

T(i) _ l + l-vTiAi) , . 

- 21 + [A,) ' ^""^^ 

with a corresponding expression for body 2, where 

Note that the overall scale of ff has no effect on the apsidal constant, to linear order. 

For a uniformly rotating body, the disturbing potential at a point x' is the centrifugal potential 



1 2 /2 47r 2 ,2 

-uo r u r 

3 15 



E^2;.(A)r2„(f7'). (B6) 



Thus we read off the rotational coefficient for body 1 



dL^'^=-J^?>^;,n(Ai) (B7) 

with the coefficient for body 2 obtained by interchange (the spherical coefficient d only contributes at second order 
in small quantities). Similarly defining — f2^Y2„^{X), we find that /|*' also satisfies Eq. (|B3ap for I = 2, but with 
the boundary condition 

+ 2/2^1) = -^^. (B8) 

3 TOl 



The rotational apsidal constant is also independent of m, and, since the overall scale of /2^(^i) is irrelevant, it is equal 
to the I = 2 tidal apsidal constant: fc^*-^-* — fcj*-^'' = This equality will not hold when non-linear corrections are 
included (it will also not hold in general relativity, when frame-dragging and other relativistic effects are included). 
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2. Energy, angular momentum and equations of motion 

23, (|X26|l . (jX27|l . (|X3T|l . and |X32)) . and making- 



Substituting these results for d^^, d^rm ^^'^ "^Im; i'^to Eqs 
use of Eq. (|A24|) , we obtain 



T = 
J* = 



L 



<jkl> 



(1) 



^Thermal + ^ 3^l'^l'^2 



fim 
r 

-(1^ 



mi"! 



2 2 i "^2.„ i i. 

-LUiUJi g-(3n a;i • n - ujI) 



+ (1^2) 



6^ + - 5"XAi • n)2 + 2AlAi • n] 



<jk> _ 



m .-7 , (^ \ mon 



(B9a) 
(B9b) 
(B9c) 

(B9d) 

(B9e) 
(B9f) 



It is simple to show that the equation of motion ljB9d|) admits the two conserved quantities 

TO2 1 m2wj 



+AlfcW^ + (1^2) 



3 r3 



-[3(Ai -n) 



J" = ^(x X v)^ + 2m2A\k'i^ujlB\{t) + (1^2), 



(BlOa) 
(BlOb) 



where B\{t) = /*(n x Ai)*(n • Ai)r ^dt, where we assume that the various quantities entering the perturbing terms 

(Ai, u)\, Ai etc.) are constant in time to the order considered. By comparing Eqs. (|B9c|) and (|B10bll . we see 
that the total constant angular momentum can be written in the form 



J* = ^(x X v)* + Jl + + ^2m2A\k^^'^ujlB\{t) + (1^2)}, 
where J\ and J| are separately constant, defined by J\ — IiCul, where the constant uj\ is given by 

"2 2. .1 "12, 



o /i5i,(i) r 

Zi -ti-'Y '•'2 



3 /i 



-(37i'a;i • n - ^) - 3m2wfBJ(t) 



(Bll) 



(B12) 



Notice that B\{t) is orthogonal to n and Ai, and vanishes if the body's spin axis is perpendicular to the orbital plane. 
Calculating the total energy E — T + W from Eqs. (|B9a(l and l|B9b|l and converting from uj\ to the constant uil, we 
obtain the final conserved energy, including tidal and rotational contributions 



5i,(l) 
2 



9 ^ 



(B13) 



Modulo constants, this is identical to E* , Eq. (|B10a|) . 
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3. Clairaut's equation and the apsidal constants 



To determine the tidal and rotational distortion effects in our binary system, it is sufficient to know the disturbing 
forces (leading to the coefficients dim) and the apsidal constants. To linear order, the apsidal constants can be obtained 
from solutions of Eq. ljB3a|l . along with Eq. ljB4p : this applies to both tidal and rotational perturbations. Because 
the scale of /; is irrelevant to the value of ki , it is useful to recast Eq. (jB3a|l into a first-order differential equation for 
77/, sometimes called Radau's equation 

077/ + 6Virji + 1) + rnim -l)-l{l + l) = 0, (B14) 

where D — Anp{a)a^ /3'm,(a) — p{a)/p{a), and where we drop the superscripts T or R. Near the origin, where 
I? — > 1, the regularity of /; requires that ry/(a) — > Z — 2 as a ^ 0. Given a density profile for a spherically symmetric 
configuration provided by a chosen equation of state, one integrates Eq. HB14p from the center to the surface, thereby 
obtaining 77; {A) , and thus fc; . 

Exact solutions of Radau's equation are known for special cases. For a homogeneous star, with p = constant, I? = 1, 
it is easy to show that rn = 1 — 2, and 

^^Homogeneous ^ . (B15) 

For a point mass, with D = except at the origin, 77/ = Z + 1, and hence, as expected, = 0. Generally, if the 

density nowhere increases outwards (i.e. if p'{a) < 0), then rii{A) satisfies the inequalities I — 2 < rii{A) < I + 1. For 
nearly homogeneous configurations, and for / = 2, Radau's equation can be rewritten in the approximate form 



da 



l + ^772(a)2 + 0(772 (a)3) 



(B16) 



Since 772(a) = in the homogeneous limit, one can take the lowest order approximation, integrate, and then, after 
some manipulation, show that 

772 « 3(1 - I/Ih) + 0[(1 - I/Inf] , (B17) 

where / is the moment of inertia of the star, and Ih is its homogeneous counterpart. Finally, for polytropic Newtonian 
stars, with equation of state p = kp^ = kp^~^^^", Kopal lists computed values of 77; and ki in Tables 2-1 and 2-2, 
for I = 1 ... 7 and n = ... 5. For 77 = 1 or F = 2, which is a common choice in numerical models of binary neutron 
stars, fc2 = 0.260 and k^ = 0.106. Kopal "2^, '2^ explores these and other general properties of Radau's equation. 

It is also useful to be able to read off values of fc2 from the external field or spacetime geometry of a given solution. 
For example, if the perturbing coefficient has the form dim — diYi*^{X), where A is a principal axis of the perturbation, 
then the external potential Eq. (|A14|I takes the form 

77? ' /l^^"'""'^ 

U=- + 2j2^j^kidiPi{cos^), (B18) 
/ 

where 7 is the angle between A and the field point. Then, if one can read off the multipole moments Qi from an 
expansion of the field in the form m/r + (5;P;(cos7)/r'+^, then one can determine the apsidal constants according 
to 

(B19) 



' 2A2'+id, ■ 

For I = 2 rotational perturbations, d2 = —uj'^/S, so fc2 — {■i/2)Q2/A^u!'^. This permits one to determine fc2 from a 
numerical solution of a rotating neutron star, for example. For Kerr black holes, Q2 = S'^/M. To lowest order in 
MiyyLu, where Mi„ is the irreducible mass, S = 4:M^^.^oj and A = 2Mi„, so that fc^^ = 3/4, precisely the value for a 
homogeneous Newtonian star. 
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